clear;clc;

FolderList      =   {'Shared','Temp','SteadyState','Equations','Graph'};
for ii=1:length(FolderList)
    addpath(genpath(FolderList{ii}));
end

%% Setup
StdVec      =   [1 1 1 -1 -1 60 -10]'*0.01/4;
InputParam  =   struct('KAPPA',0.00/4,'XI',1/82.5/4,'Pir_SS',0.02/4,...
                                'Flag_dE',0,'Flag_Price',0,...
                                'B_Total_SS',0.91,'b_lb',-0.29,...
                                'TrProb_H',0.98,'Prob_H',0.37,...
                                'TrProb_ext',0.92,'Prob_ext',0.33,...
                                'Fin_ExtBond',0.1930,'Fin_AdjCost',0.8,...
                                'Cons_FracH',0.60,...
                                'RHO_M_dom',0.68,'RHO_M_ext',0.81,'RHO_Y_H',0.5,'RHO_p_F',0.90,...
                                'Taylor_Pi',1.1,'Taylor_ir',0.87,...
                                'Cons_ElasTN',6.19,'Cons_ElasHF',6.19);
PP          =   Setup_PP(InputParam);


%% Steady State
[SS,PP]     =   SteadyState_Solver(PP);

%% Construct the Model Structure
MODEL       =   Setup_MODEL(PP,SS);

% Check the Steady State in Equation Blocks
Res_NK      =   Equ_NK(PP,SS,0);
Res_Portfolio=  Equ_Portfolio(PP,SS,MODEL,0);
Res_HH      =   Equ_HH(PP,SS,MODEL,0);

if max(abs([Res_NK;Res_Portfolio;Res_HH]))<1e-5
    fprintf('Steady State is Correct!\n');
else
    error('Steaty State is not Correct, check it  ...');
end


%% Compute the Jacobian
EquJac          =   struct();
EquJac.HH       =   Equ_HH(PP,SS,MODEL,1,[],[],1);
EquJac.HH_Aux   =   Equ_HH_Aux(PP,SS,MODEL,1,[],[],1);
EquJac.NK       =   Equ_NK(PP,SS,1);
EquJac.Portfolio=   Equ_Portfolio(PP,SS,MODEL,1,[],[],1);

%% IRF
% Flexible Exchange Rate
PP.Taylor_dE    =   0;
PP.Fin_AdjCost  =   0.1;

EquJac.NK       =   Equ_NK(PP,SS,1);
EquJac.Portfolio=   Equ_Portfolio(PP,SS,MODEL,1,[],[],1);
fir_ord         =   JacAssemble(MODEL,EquJac);
tic;
[gx,hx,gxhx_ExitFlag,gxhx_Info]=gx_hx(fir_ord,1);
fprintf(['Perturbation is Solved, Elapse Time=',num2str(toc,'%.2f'),'s\n']);
if gxhx_ExitFlag~=1
    warning('Perturbation Solution is problematic');
    SOLUTION    =   [];
    return
else
    solution    =   struct('gx',gx,'hx',hx);
end
IRF             =   IRF_1order(MODEL,solution,StdVec,[],4*10);
IRF             =   TempFun_AugIRF(PP,SS,MODEL,EquJac,IRF);

%% Comparative statics 
PP.Fin_AdjCost          =   0.01;
[IRF_LL, SOLUTION_LL]   =   TempFun_NewIRF(PP,SS,MODEL,EquJac,StdVec);
PP.Fin_AdjCost          =   0.1;
[IRF_L, SOLUTION_L]     =   TempFun_NewIRF(PP,SS,MODEL,EquJac,StdVec);
PP.Fin_AdjCost          =   0.5;
[IRF_H, SOLUTION_H]     =   TempFun_NewIRF(PP,SS,MODEL,EquJac,StdVec);


SubFun_CollectResults(SOLUTION_L)


save('SOLUTION.mat', 'SOLUTION_LL', 'SOLUTION_L', 'SOLUTION_H');
save('IRF.mat', 'IRF_LL', 'IRF_L', 'IRF_H');

%% Key Results 
TempPlot    =	@(xx)plot(xx,'linewidth',3);
Folder      =   'TableGraph/';
% External monetary shock:
% Aggregate C
figure; 
TempPlot(IRF_LL.Eps_Eta_M_ext.C*100); hold on
TempPlot(IRF_L.Eps_Eta_M_ext.C*100); hold on
TempPlot(IRF_H.Eps_Eta_M_ext.C*100);
legend("adj cost=0.01", 'adj cost=0.1', 'adj cost=0.5');
print('-dpng','-r1000', [Folder,'ExtMon_AggC']);
% C-Diff: ext - dom
figure; 
TempPlot(diff(IRF_LL.Eps_Eta_M_ext.AggStat(1:2, :))*100); hold on
TempPlot(diff(IRF_L.Eps_Eta_M_ext.AggStat(1:2, :))*100); hold on
TempPlot(diff(IRF_H.Eps_Eta_M_ext.AggStat(1:2, :))*100); 
legend("adj cost=0.01", 'adj cost=0.1', 'adj cost=0.5');
print('-dpng','-r1000', [Folder,'ExtMon_CDiff']);
% UIP dev
figure; 
TempPlot(IRF_LL.Eps_Eta_M_ext.DevUIP*400); hold on
TempPlot(IRF_L.Eps_Eta_M_ext.DevUIP*400); hold on
TempPlot(IRF_H.Eps_Eta_M_ext.DevUIP*400);
legend("adj cost=0.01", 'adj cost=0.1', 'adj cost=0.5');
print('-dpng','-r1000', [Folder,'ExtMon_UIPDev']);
% 
% figure; 
% TempPlot(IRF_LL.Eps_Eta_M_ext.ir_ext*400); hold on
% TempPlot(IRF_L.Eps_Eta_M_ext.ir_ext*400); hold on
% TempPlot(IRF_H.Eps_Eta_M_ext.ir_ext*400);
% legend("adj cost=0.01", 'adj cost=0.1', 'adj cost=0.5');

% Domestic monetary shock: UIP dev
figure; 
TempPlot(IRF_LL.Eps_Eta_M_dom.DevUIP*400); hold on
TempPlot(IRF_L.Eps_Eta_M_dom.DevUIP*400); hold on
TempPlot(IRF_H.Eps_Eta_M_dom.DevUIP*400);
legend("adj cost=0.01", 'adj cost=0.1', 'adj cost=0.5');
print('-dpng','-r1000', [Folder,'DomMon_UIPDev']);


% Key results 
TempList    =   {SOLUTION_LL, SOLUTION_L, SOLUTION_H};
LabelList   =   {'LL','L','H'};

TempNum     =   length(TempList);
Tables      =   cell(TempNum,1);

parfor ii=1:TempNum
    vv          =   TempList{ii};
    Tables{ii}  =   SubFun_CollectResults(vv);
end

FileName    =   'TableGraph\\ModelExtension_AdjCost.xlsx';
ShockList   =   {'Eps_Eta_M_dom', 'Eps_Eta_Y_H', 'Eps_Eta_M_ext'};
for ii=1:length(ShockList)
    SH          =   ShockList{ii};
    TempData    =   [];
    for jj=1:TempNum
        TempData    =   [TempData,Tables{jj}.(SH)];
    end
    TempTable   =   array2table(TempData,'VariableNames',LabelList,'RowNames',Tables{1}.Row);
    writetable(TempTable,FileName,'WriteRowNames',true,'Sheet',SH);
end
